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ABSTRACT 

Spherical solar dynamo simulations are performed. Self-consistent, fully compressible magnetohy- 
drodynamic system with a stably stratified layer below the convection zone are numerically solved with 
a newly developed simulation code based on the Yin- Yang grid. The effects of penetrative convection 
are studied by comparing two models with and without the stable layer. A solar-like differential rota- 
tion profile is established when the penetrative convection is taken into account without assuming any 
forcing. A large-scale magnetic field is also spontaneously organized in the underlying stable layer. 
The embedded field has a dipole symmetry about the equator and it shows polarity reversals in time. 
Subject headings: convection-magnetohydrodynamics (MHD) - Sun; dynamo - Sun; interior 



1. INTRODUCTION 

Solar interior is a frontier of the solar physics. A 
grand challenge is a construction of self-consistent theory 
that explains the observed large-scale spatial structures 
of the fields and their dynamical change in time with 
the pseudo-periodicity of 22 years. Two basic large-scale 
structures that should be explained are the azimuthal av- 
erage of the azimuthal flow, w^, and the azimuthal aver- 
age of the azimuthal magnetic field, -B^. See lOssendriiven 
[2003; Miesch 2012 for reviews of the solar interior. 

To reproduce the large-scale structures and dynamics, 
magnetohydrodynamics (MHD) simulations have been 
perfo rmed both in the global (spherical shell) geometr y 
(e.g., iGilmanlHOSl 'Glatzmaier|[l98l; |Brun et al."2003 ) 
and i n the local Cart esian ge ome t ry (e.g.,|N ordlund et ajj 
Il992t iBrandenbur'g et al.l ll"99l ICattaneo et all I2003D . 
The exponential growth of the computer performance 
has yielded br eakthroughs to the sola r dynamo simula- 
tion research. iBrowning et al.l ()2006l) showed in spher- 
ical shell dynamo simulations that strong axisymmetric 
toroidal magnetic fields can be formed in stably strat- 
ified layer below the convection zone. A latitudinal 
thermal forcing was, however, necessary for keeping the 
solar-like different ial rotation profile in their simulation. 
IB run et al.l (|2011l ) performed hydrodynamic simulations, 
without magnetic fields, in which the solar-like convec- 
tion profile with a tachocline can be formed without any 
forcing. An anelastic model was adopted in these two 
simu lations. 

Gh izaru et all (|2010D and lRacine et al.l (|201lD showed 
that the large-scale magnetic field in the stable layer ex- 
hibits polarity reversals when the temporal integration 
of the simulation is calculated for long enough. Their 
simulation is also based on an anelastic model that is 
popularly used in the global circulation models of the 
earth's atmosphere with a cooling term to force the sys- 
tem toward the ambient state. 

In this Letter, we report the first fully compressible 
spherical solar dynamo simulation with a stably strati- 
fied layer below the convection zone. Formations of the 
key profile of the solar interior, i.e., the solar-like Vcj, & 
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i?0, and the spontaneous polarity reversals are simulated, 
without assuming any forcing in the fundamental equa- 
tions. To elucidate the effects of the penetrative convec- 
tion, two simulations with and without the stable layer 
below the convection zone are compared. 

Another purpose of this Letter is to report a develop- 
ment of new program code for the solar dynamo simu- 
lation. Most simulation models of the global solar dy- 
namo are spectral-base d type, using th e spherical har- 
monics expansion (e.g.. lBrun et al.ll2004[ ). The spherical 
harmonics expansion method is, however, believed to be 
confronted with the parallel scaling difficulty when tens 
of thousands, or more, processor cores are used. A differ- 
ent approach to massively parallel solar dynamo model 
is strongly required for the present peta- or coming exa- 
scale era. We have developed a global solar dynamo sim- 
ulation code based on the grid point-based approach. 

The spherical geometry imposes difficulties in the 
design of the spatial grid points to sustain high nu- 
merical efficiency, accuracy, and parallel scalability. 
We have proposed an overs et grid method ap proach 
to the spherical geometry (jKagevama fc Sata [2004) . 
The grid system, Yin- Ya ng grid, is applied to geo- 
dyna mo (e.g., Kageva ma et al.l 120081: iMivaeoshi et al.l 
20101). mantle convection (e.g., iKamevama et al.l 120081: 
Tacklevll2008D. supern ova explosions (e.g., iMiiller et al] 
2012tlLentz et al.ll2012D . and other astro- and geophysical 
simulations. The parallel scaling property of the spheri- 
cal MHD simulation on the Yin- Yang grid is promising. 
It attained 46% (15.2 TFlops) of the peak performance 
of 4096 cores of the Earth Simulator supercomputer for 
the geodynamo simulations (Gordon Bell Award in Su- 
percomputing 2004). Our new solar dynamo code is de- 
veloped based on this Yin- Yang geodynamo code. This 
Letter is our first report on the results obtained by this 
Yin- Yang solar dynamo code. 

2. NUMERICAL SETTINGS 

We numerically solve an MHD dynamo convection sys- 
tem in a spherical shell domain defined by (0.6i? < r < 
R), (0 < < tt), and (— tt < < tt), where r, 9, and 
are radius, colatitude, and longitude, respectively. Our 
model has two layers: upper convective layer of thickness 
0.3R in the range of {0.7 R < r < R), and stably stratified 
lower layer of thickness O.IR in {0.6R < r < 0.7R). 
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Fig. 1. — The radial velocity profile Vr{r,8) when t = 330tc, in the MoUweide projection, on spherical surfaces at different depths for 
two models A and B. The panels (a)-(c) correspond to the depths r = 0.95-R, 0.85-R and 0.72R for Model A, and the panels (d)-(f) are 
those for Model B. The orange and blue tones depict upflow and downflow velocities. 



The fundamental equations are the fully compressible 
MHD equations in the rotating frame of reference with 
a constant angular velocity CI which is parallel to the 
coordinate axis (6 = 0): 



dt 
df 

dt 



= —V • 



/, (1) 

(vf) -Vp + jxB + pg + 2pvxn + Q,{2) 
\7p - -fp\/ • t; + (7 - 1)$ , (3) 

(4) 



dp 

dA . 

— =v X B -T]j 

with g = —go/r^Br, B = VxA, j — VxB, and 
$ = V • (kVT) -I- 2^8'^ + r]j^. The mass density p, pres- 
sure p, mass flux density / = pv, magnetic field's vector 
potential A are the basic variables. We assume an ideal 
gas law p = (7— l)/5e with 7 = 5/3, where e is the internal 
energy. The viscosity, magnetic diffusivity, and thermal 
conductivity are represented by ly, i], and k respectively. 
The viscous term Q is written by Qi = dj {2vSij) with 
the strain tensor 



dxj 



dxi 



r%v 



(5) 



In the program code, the tensor components are written 
in the spherical polar coordinate. 

The initial condition is a hydrostatic equilibrium which 
is described by a piecewise polytropic distribution with 
the polytropic index to, 

(6) 



dr Cv (7 
iKapvla et aP [2010l) . 



l)(m+l) 



(e.g., iKapvla et al.l 120100 . We choose to = 1 and 3 
for the upper convection layer and the lower stable 
layer, respectively. The thermal conductivity is deter- 
mined by requiring a constant luminosity L, defined by 
L = —AiTKr^dT/dr, throughout the domain. 

Non-dimensional quantities are defined by setting R = 
ffo = Po = 1 where po is the initial density at r = 0.6R. 
We normalize length, time, velocity, density, and mag- 
netic field in units of R, ^/RFJgo, ^^go/R, po and 



^/gopo/R. We define the Prandtl, magnetic Prandtl, and 
Rayleigh numbers by 



Pr=-, Pm==-, Ra 



X 



n 






Cp dr 



(7) 



where x = /"^m/icpPm) is the thermal diffusivity at the 
mid-convection zone (r = rm), and d = 0.3-R is the depth 
of the convection zone. The stratification level is con- 
trolled by the normalized pressure scale height at the 
surface, ^0 = Cv{j — l)Ts/{go/R), where Tg is the tem- 
perature at r = i?. In this work, we use ^0 = 0-3, yielding 
a small density contrast about 3. 

The relative importance of rotation in the convection 
is measured by the Coriolis number defined by Co = 
2f7o<^/i'rms, where Wi-ms = [{vg + vf.)]^^'^ is the mean ve- 
locity. The angular brackets denote the time and volume 
average in the upper convective layer at the saturated 
state. The convective turn-over time and the equiparti- 
tion strength of magnetic field are defined by Tc = d/vrms 

The stress-free boundary condition for the velocity is 
imposed on the two spherical boundaries. We assume 
the perfect conductor boundary condition for the mag- 
netic field {Ar — Ag — A^ = 0) on the inner sur- 
face, and the radial field condition {A^ = 0, dAg/dr — 
—Ag/r,dA^/dr — —A^/r) on the outer surface. A con- 
stant energy flux is imposed on the inner boundary. The 
temperature is fixed to be Tg on the outer boundary. 

The eqs. ([T])-(l4]) are discretized by the second-order 
central difference on the Yin- Yang grid. Each of the 
two congruent grids, Yin-grid and Yang grid, covers a 
partial spherical shell region defined as (7r/4 < < 
37r/4, — 37r/4 < 4> < 37r/4). They are combined in a com- 
plemental way to cover a whole spherical shell. Physical 
quantities on the horizontal boarders of Yin- or Yang- 
grid are set by mutual interpolations. For the time inte- 
gration, the standard fourth-order Runge-Kutta method 
is used. Since the Yin- Yang grid is free from the coordi- 
nate singularity and the grid concentration around there, 
we can avoid the severe time-step constraint due to the 
CFL condition. See Kagevama fc Sato..2004 for more de- 
tail on the Yin- Yang grid method. The computation is 
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performed in parallel using MPI (Message Passing Inter- 
face) . 

Non-dimensional parameters Pr = 0.2, Pm — 4.0, and 
Ra = 1.2 X 10^, and constant angular velocity of i7o — 0.4 
are adopted in all the calculations reported here. The 
total grid size for the run with the upper convection layer 
and the lower stable layer (Model A) is 121 (in r) x402 
(in Q) X 402 (in (/)) x2 (Yin & Yang). A model without 
the stable layer (Model B) is also studied, in the domain 
{0.7 R < r < R), with the same physical parameters and 
the same grid spacings (91 x 402 x 402 x 2). A random 
temperature perturbation and weak magnetic field are 
seeded in the convection zone when the calculation starts. 

3. NUMERICAL RESULTS 

3.1. Convection and Mean Flows 

After the convective motion sets in, it reaches a nonlin- 
ear saturation state at around t = IOtc. The saturation 
levels of the convection kinetic energy for Model A and B 
are almost the same. The mean velocity is Wrms — 0.03 
which yields Bcq = 0.02, Co = 8.0 and r^ = 10.0 for 
both models. We have run the simulations till 500tc and 
then compare physical properties of convections, mean 
flows and magnetic dynamos between two models. The 
flow velocity and magnetic flux density are normalized 
by ?^rms and Beq respectively in the following. 

Figure 1 shows, in the Mollweide projection, the ra- 
dial velocity profile Vr when t = 330tc on spherical 
surfaces at different depths for two models. The pan- 
els (a)-(c) correspond to the depths r — 0.95i?, 0.85i? 
and 0.72i? for Model A, and the panels (d)-(f) are those 
for Model B. The orange and blue tones depict upflow 
and downflow velocities. At the upper (r = 0.95i?) and 
mid (r = 0.85i?) convection zones, the convective mo- 
tion is characterized by upflow dominant cells surrounded 
by networks of narrow downflow lanes for both models. 
The higher the latitude, the smaller the convective cell 
prevails. Elongated columnar convective cells aligned 
with the rotation axis appear near the equator. These 
are the typical f eatures obs erved in rotating stratified 
convection (e.g. , ISoruit et a l. 1990; Miesch ct al. 2000; 
iBrummell et al.l |2002|) . In panel (c) , we find that the 
downflow lanes persist the plume-like coherent struc- 
ture even just above the bottom of the unstable layer 
(r = 0.72i?). The downflow plumes then penetrate into 
the underlying stable layer. As a result, mean zonal and 
meridional flows are driven by the Reynolds stress in the 
stable layer. 

In Figures 2(a) and (b), time-averaged mean angular 
velocity, defined by Vt{r^6) = u^/(rsin0) -|- i7o, is shown 
for the models A and B, respectively. Here the overbar 
denotes the time and azimuthal average. The time aver- 
age spans in the range of 300tc < t < 400rc. While the 
overall structures -especially the equatorial acceleration- 
are common in the two models A and B, some differences 
between them are noticeable. Model B exhibits a cylin- 
drical rotation profile which has iso-rotation contours 
aligned with the rotation axis. The system is dominated 
by th e Taylor-Proudman balance in this case (jPedlosky I 
|1987() . In contrast to that, a differential rotation with 
conical iso-rotation surfaces at mid-latitude is established 
in Model A. More importantly, a strong radial gradient of 
the angular velocity is developed in the stably stratified 
layer around latitudes ±40°. This structure is reminis- 



cent of the solar tachocline despite the rad ial shear layer 
is br oad compared to the observed one (jHughes et al.l 
|2007( ). The rotation profile of Model A is reasonably 
simil ar with that of the Su n deduced from helioseismol- 
ogy ([Thompson et al.1l2003[ l. This confirms that the sta- 
ble layer underlying the convective envelope has a role 
in maintaining the solar-like differen tial rotation as re - 
ported in the hydrodynamic model of lBrun et al.l ()2011fl . 
Shown in Figures 2(c) and (d) are time-averaged mean 
meridional flows for the models A and B. The color con- 
tour depicts the absolute value of the meridional flow 
velocity with a maximum ~ O.lwrms- Overplotted are 
streaklines with a length proportional to the flow speed. 
The circulation flow is primarily counter-clockwise in the 
bulk of the convection zone in the northern hemisphere, 
that is, the poleward in the upper convection zone and 
the equatorward in the bottom convection zone in both 
models. There is however a clear difference in the cir- 
culation pattern between the two models. While a large 
single-cell is formed in Model B, Model A shows a double- 
cell pattern with a strong inward/outward flow at the 
low/mid latitudes as suggested by the hel ioseismic inver- 
sion (see lMitra-Kraev fc Thompsonll2007l ) . An intriguing 
finding is that the equatorward component penetrates 
into the underlying stable layer when the strong radial 
gradient of the angular velocity resides. The penetrative 
transport of magnetic flux by the meridional flow might 
play a role in magnetic dynamo in our model. 

3.2. Structure and Evolution of Magnetic Fields 

The magnetic energy is amplified by the dynamo action 
and is saturated at a level of about 40% of the convective 
kinetic energy for both models after t ~ 20tc. A snap- 
shot of the azimuthal component of the magnetic field 
i?0 at i = 330tc is presented in Figure 3 on a spherical 
surface at the depths (a) r = 0.62i? and (b) r = 0.85i? 
for Model A, and (c) r = 0.85i? for Model B. The orange 
and blue tones depict positive and negative values of the 
B^ component. The magnetic field lines at the time and 
position corresponding to those in the panels (a)-(c) are 
visualized in Figures 3(d)-(f), respectively. 

The convective envelope is dominated by disordered 
tangled magnetic field lines with a myriad of localized 
small-scale structures in both models. These incoher- 
ent magnetic fields are strongly infiuenced by vigorous 
convective motions and thus are highly intermittent. 
The horizontal converging flows sweep magnetic flelds 
into downflow lanes and intensify them locally to the 
super-equipartition strength a s was observed in convec- 
tive dynamo simulations fe.g.. 'Bran denburg et al.lll996l : 
Cattanco ct al. 2003; Brun ct al. 200^^ 

In the underlying stable layer of Model A, a strong 
large-scale azimuthal component is built up and stored 
for long time intervals. This well-organized magnetic 
component is roughly antisymmetric about the equa- 
torial plane and reaches a maximum strength of 2i?oq. 
The large-scale component grows in the region where the 
strong radial shear resides (see Figure 2(a)). This would 
be an important evidence of a connection between the 
deep-seated coherent fields and the strong radial shear 
layer. These results suggest that the stably stratified 
layer is a key component to organize large-scal e magnetic 
fields and support numer i cal st udies of iBrowning et aLl 
(|200^ and lGhizarueTall (l2010l) . 
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Fig. 2. — Time-averaged mean angular velocity Q(r,8) [panels (a) and (b) for Models A and B], and mean meridional flow [panels (c) 
and (d) for Models A and B], where Q{r,8) = ^^/(rsinS) + Qq. The white solid curves in panel(a) and (c) denote the interface between 
the convection and stable layers. 



One of the most interesting findings in our simulation 
is that the deep-seated magnetic fields in the stable layer 
show polarity reversals. Figure 4 gives an azimuthally- 
averaged magnetic field as a function of time and lat- 
itude._ The panels (a), (b) and (c) represent Br, B_e, 
and B^ at r = 0.65R for Model A. The large-scale B^ 
with antisymmetric parity persists over a relatively long 
period despite strong stochastic disturbances due to pen- 
etrative convective motions. The B^ component changes 
the sign for at lest three times, at about t = IOOtc, 
t — 200tc, and t = SSOtc [panel (c)]. 

As for the poloidal component of the magnetic field, 
it shows_ the dipole dominance. During 30-100tc when 
strong B^ component with positive polarity dominates 
in the northern hemisphere, negative Br and positive Bg 
are observed [panels (a) and (b)]. While the Bg compo- 
nent has the same sign in both the hemispheres, the Br 
component Jras opposite polarity in the two hemispheres. 
When the B^ reversal takes place, the other two compo- 
nents Br and Bg also change the sign. See, for example, 
at about t = IOOtc in Fig. Ill 

The overall features indicate that the large-scale mag- 
netic field has a dipole-like component rooted in the sta- 
bly stratified layer and reverses its dipole polarity, ac- 
companied by the polarity reversal of the B^ component 
as shown in Figure 5 schematically. 



4. SUMMARY 

We reported, in this Letter, our first results of solar 
dynamo simulation based on the Yin- Yang grid with the 
fully compressible MHD model. To investigate influences 
of the stably stratified layer below the convection zone, 
two simulation models with and without the stable layer 
were compared. 

It is confirmed that the stable layer has substantial in- 
fluence on the convection and the magnetic field. The 
differential rotation profile inferred from the helioseis- 
mology is simulated when the stable layer is taken into 
account. The conical iso-rotation surfaces at the mid 
latitudes is established in the convective envelope. The 
strong radial shear layer, which is reminiscent of the solar 
tachocline, is formed just beneath the convection zone as 
a result of the penetrative convection. 

The stably stratified layer would be a key component 
not only for maintaining the solar-like differential rota- 
tion but also for organizing large-scale magnetic fields. In 
our model, the strong large-scale azimuthal component 
with antisymmetric parity about the equator is built up 
in the underlying stable layer, whereas the convective 
envelope is dominated by highly intermittent disordered 
magnetic fields. An intriguing finding is that the large- 
scale magnetic fields embedded in the stable layer un- 
dergo polarity reversals. Although a full description of 
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Fig. 3. — A snapshot of the azimuthal component of the magnetic field B,p when t = 330tc on a spherical surface at the depths (a) 
r = 0.62R and (b) r = 0.85-R for Model A, and (c) r = 0.85-R for Model B. The orange and blue tones depict positive and negative values 
of the B^ component. The magnetic field lines at the time and position corresponding to those in the panels (a)— (c) are visualized in the 
panels (d)-(f), respectively 
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_FlG_. 4. — An azimuthally-averaged magnetic field as a function of time and latitude. The top, middle and bottom panels correspond to 
Br, Bg, and Bj, at the depths r = 0.65-R in Model A. The red and blue tones depict positive and negative values of each component. 
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Fig. 5. — Schematic picture of the large-scale magnetic field 
deep-seated in the convectively stable layer. The blue and red 
curves demonstrate the azimuthal and poloidal magnetic compo- 
nents. The large-scale magnetic field has a dipole-like component 
and reverses its dipole polarity, accompanied by the polarity rever- 
sal of the B^ component. 

the dynamo mechanism is beyond the scope of this Let- 
ter, here we mention that the strong radial shear layer 



developed just beneath the convection zone is closely re- 
lated to the dynamo process. 

All the dynamo simulations reported here have used 
a relatively weak stratification with the density contrast 
of about 3 (see §2). The strong stratification in the real 
Sun may influence on the physical propert ies of convec- 
tions , mean flows and magnetic dynamo ([Kapvla et al.l 
|2013| ). It would be interesting that the three key fea- 
tures, solar-like v^, B^, and the polarity reversals are 
self-consistently reproduced, without assuming any forc- 
ing, even in the modest density stratification. Higher 
resolution simulations with a more realistic density strat- 
ification will facilitate our understanding of the physics 
of the solar convection and the solar dynamo, that is our 
next step with the Yin- Yang solar dynamo simulation 
code. 
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